home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / pascal / tpl60n19.zip / TESTPRGS.ZIP / UNIT2.PAS < prev    next >
Pascal/Delphi Source File  |  1992-02-27  |  25KB  |  822 lines
  1. {$a+,n-,x-,s-,i-,r-,b-,v-}
  2.  
  3. unit Unit2;
  4. interface
  5.    uses mainvars;
  6.    procedure mile70170;
  7. implementation
  8.    procedure mile70170;
  9.    begin
  10.  
  11. {=============================================}
  12.    Milestone := 70;
  13. {=============================================}
  14.    writeln;
  15.    writeln ('Running test of square root(x).');
  16.    TestCondition (Failure, (Zero = sqrt (Zero))
  17.          and (- Zero = sqrt (- Zero))
  18.          and (One = sqrt (One)), ' Square root of 0.0, -0.0 or 1.0 wrong  '
  19.          );
  20.    MinSqrtError := Zero;
  21.    MaxSqrtError := Zero;
  22.    J := 0;
  23.    X := Radix;
  24.    OneUlp := U2;
  25.    SqrtXMinX (SeriousDefect);
  26.    X := BInverse;
  27.    OneUlp := BInverse * U1;
  28.    SqrtXMinX (SeriousDefect);
  29.    X := U1;
  30.    OneUlp := U1 * U1;
  31.    SqrtXMinX (SeriousDefect);
  32.    if J <> 0 then
  33.       begin
  34.       NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
  35.       Pause;
  36.       end;
  37.    writeln ('Testing if sqrt(X * X) = X for ', NoTrials, ' integers X.');
  38.    J := 0;
  39.    X := Two;
  40.    Y := Radix;
  41.    if (Radix <> One) then
  42.       repeat
  43.          X := Y;
  44.          Y := Radix * Y;
  45.       until (Y - X >= NoTrials);
  46.    OneUlp := X * U2;
  47.    I := 1;
  48.    Continue := true;
  49.    while (I <= NoTrials) and Continue do (* dgh: 10 --> NoTrials *)
  50.       begin
  51.       X := X + One;
  52.       SqrtXMinX (Defect);
  53.       if J > 0 then
  54.          begin
  55.          Continue := false;
  56.          NoErrors [Defect] := NoErrors [Defect] + 1;
  57.          end;
  58.       I := I + 1;
  59.       end;
  60.    writeln ('Test for Sqrt Monotonicity.');
  61.    I := - 1;
  62.    X := BMinusU2;
  63.    Y := Radix;
  64.    Z := Radix + Radix * U2;
  65.    NotMonot := false;
  66.    Monot := false;
  67.    while not (NotMonot or Monot) do
  68.       begin
  69.       I := I + 1;
  70.       X := sqrt (X);
  71.       Q := sqrt (Y);
  72.       Z := sqrt (Z);
  73.       if (X > Q) or (Q > Z) then
  74.          NotMonot := true
  75.       else
  76.          begin
  77.          Q := Int (Q + Half);
  78.          if (I > 0) or (Radix = Q * Q) then
  79.             Monot := true
  80.          else if I > 0 then
  81.             begin
  82.             if I > 1 then
  83.                Monot := true
  84.             else
  85.                begin
  86.                Y := Y * BInverse;
  87.                X := Y - U1;
  88.                Z := Y + U1;
  89.                end
  90.             end
  91.          else
  92.             begin
  93.             Y := Q;
  94.             X := Y - U2;
  95.             Z := Y + U2;
  96.             end
  97.          end
  98.       end;
  99.    if Monot then
  100.       writeln ('Sqrt has passed a test for Monotonicity.')
  101.    else
  102.       begin
  103.       NoErrors [Defect] := NoErrors [Defect] + 1;
  104.       writeln('DEFECT:  Sqrt(X) is non-monotonic for X near ', Y);
  105.       end;
  106. {=============================================}
  107.    Milestone := 80;
  108. {=============================================}
  109.    MinSqrtError := MinSqrtError + Half;
  110.    MaxSqrtError := MaxSqrtError - Half;
  111.    Y := (sqrt (One + U2) - One) / U2;
  112.    SqrtError := (Y - One) + U2 / Eight;
  113.    if SqrtError > MaxSqrtError then
  114.       MaxSqrtError := SqrtError;
  115.    SqrtError := Y + U2 / Eight;
  116.    if SqrtError < MinSqrtError then
  117.       MinSqrtError := SqrtError;
  118.    Y := ((sqrt (F9) - U2) - (One - U2)) / U1;
  119.    SqrtError := Y + U1 / Eight;
  120.    if SqrtError > MaxSqrtError then
  121.       MaxSqrtError := SqrtError;
  122.    SqrtError := (Y + One) + U1 / Eight;
  123.    if SqrtError < MinSqrtError then
  124.       MinSqrtError := SqrtError;
  125.    OneUlp := U2;
  126.    X := OneUlp;
  127.    for Index := 1 to 3 do
  128.       begin
  129.       Y := sqrt ((X + U1 + X) + F9);
  130.       Y := ((Y - U2) - ((One - U2) + X)) / OneUlp;
  131.       Z := ((U1 - X) + F9) * Half * X * X / OneUlp;
  132.       SqrtError := (Y + Half) + Z;
  133.       if SqrtError < MinSqrtError then
  134.          MinSqrtError := SqrtError;
  135.       SqrtError := (Y - Half) + Z;
  136.       if SqrtError > MaxSqrtError then
  137.          MaxSqrtError := SqrtError;
  138.       if ((Index = 1) or (Index = 3)) then
  139.          X := OneUlp * Sign (X) * Int (Eight / (Nine * sqrt (OneUlp)))
  140.       else
  141.          begin
  142.          OneUlp := U1;
  143.          X := - OneUlp;
  144.          end;
  145.       end;
  146. {=============================================}
  147.    Milestone := 85;
  148. {=============================================}
  149.    SquareRootWrong := false;
  150.    AnomolousArithmetic := false;
  151.    RSqrt := Other; (* ~dgh *)
  152.    if Radix <> One then
  153.       begin
  154.       writeln ('Testing whether sqrt is rounded or chopped: ');
  155.       D := Int (Half + Power (Radix, One + Precision - Int (Precision)))
  156.          ;
  157.    { ... = Radix^(1 + fract) if Precision = integer + fract. }
  158.       X := D / Radix;
  159.       Y := D / A1;
  160.       if (X <> Int (X)) or (Y <> Int (Y)) then
  161.          begin
  162.          AnomolousArithmetic := true;
  163.          end
  164.       else
  165.          begin
  166.          X := Zero;
  167.          Z2 := X;
  168.          Y := One;
  169.          Y2 := Y;
  170.          Z1 := Radix - One;
  171.          FourD := Four * D;
  172.          repeat
  173.             if Y2 > Z2 then
  174.                begin
  175.                Q := Radix;
  176.                Y1 := Y;
  177.                repeat
  178.                   X1 := abs (Q + Int (Half - Q / Y1) * Y1);
  179.                   Q := Y1;
  180.                   Y1 := X1;
  181.                until X1 <= Zero;
  182.                if Q <= One then
  183.                   begin
  184.                   Z2 := Y2;
  185.                   Z := Y;
  186.                   end;
  187.                end;
  188.             Y := Y + Two;
  189.             X := X + Eight;
  190.             Y2 := Y2 + X;
  191.             if Y2 >= FourD then
  192.                Y2 := Y2 - FourD;
  193.          until Y >= D;
  194.          X8 := FourD - Z2;
  195.          Q := (X8 + Z * Z) / FourD;
  196.          X8 := X8 / Eight;
  197.          if Q <> Int (Q) then
  198.             AnomolousArithmetic := true
  199.          else
  200.             begin
  201.             Break := false;
  202.             repeat
  203.                X := Z1 * Z;
  204.                X := X - Int (X / Radix) * Radix;
  205.                if X = One then
  206.                   Break := true
  207.                else
  208.                   Z1 := Z1 - One;
  209.             until Break or (Z1 <= 0);
  210.             if (Z1 <= 0) and (not Break) then
  211.                AnomolousArithmetic := true
  212.             else
  213.                begin
  214.                if Z1 > RadixD2 then
  215.                   Z1 := Z1 - Radix;
  216.                repeat
  217.                   NewD;
  218.                until U2 * D >= F9;
  219.                if D * Radix - D <> W - D then
  220.                   AnomolousArithmetic := true
  221.                else
  222.                   begin
  223.                   Z2 := D;
  224.                   I := 0;
  225.                   Y := D + (One + Z) * Half;
  226.                   X := D + Z + Q;
  227.                   SubRout3750;
  228.                   Y := D + (One - Z) * Half + D;
  229.                   X := D - Z + D;
  230.                   X := X + Q + X;
  231.                   SubRout3750;
  232.                   NewD;
  233.                   if D - Z2 <> W - Z2 then
  234.                      AnomolousArithmetic := true
  235.                   else
  236.                      begin
  237.                      Y := (D - Z2) + (Z2 + (One - Z) * Half);
  238.                      X := (D - Z2) + (Z2 - Z + Q);
  239.                      SubRout3750;
  240.                      Y := (One + Z) * Half;
  241.                      X := Q;
  242.                      SubRout3750;
  243.                      if I = 0 then
  244.                         AnomolousArithmetic := true;
  245.                      end
  246.                   end
  247.                end
  248.             end
  249.          end;
  250.       if (I = 0) or AnomolousArithmetic then
  251.          begin
  252.          NoErrors [Failure] := NoErrors [Failure] + 1;
  253.          writeln ('FAILURE:  Anomolous arithmetic with ',
  254.             'integer < Radix^Precision = ');
  255.          writeln (W, '  fails test whether sqrt rounds or chops.');
  256.          SquareRootWrong := true;
  257.          end
  258.       end;
  259.    if not AnomolousArithmetic then
  260.       begin
  261.       if not ((MinSqrtError < 0) or (MaxSqrtError > 0)) then
  262.          begin
  263.          RSqrt := Rounded;
  264.          writeln ('Square root appears to be correctly rounded.');
  265.          end
  266.       else if (MaxSqrtError + U2 > U2 - Half) or (MinSqrtError > Half)
  267.             or (MinSqrtError + Radix < Half) then
  268.          SquareRootWrong := true
  269.       else
  270.          begin
  271.          RSqrt := Chopped;
  272.          writeln ('Square root appears to be chopped.');
  273.          end;
  274.       end